package scu.maqiang.cte;

import java.io.FileWriter;
import java.io.IOException;
import java.util.Arrays;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2Q41;
import scu.maqiang.fes.FES2Q42;
import scu.maqiang.mesh.Mesh2Q4;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.mesh.Tecplot;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class CylinderResultQ4A_SC1 {

	public static void main(String[] args) {
		long begin = System.nanoTime();
		int lineNum = 100;
		double dis = 1.0 / lineNum;
		double[] lineX = new double[lineNum + 1];
		double[] lineY = lineX.clone();
		Arrays.setAll(lineX, i -> i * dis);

		double kxx = 1.0, density = 1.0, specificHeat = 1.0;
		double E = 2.0, Nu = 0.3, alpha = 0.05;
		double beta = E / (1 - 2 * Nu) * alpha;

		double theta = 0.5, endTime = 1.0e1, dt = endTime / 1000;

		String folder = "CylinderResultQ4A\\";
		
		String heatOutputFile = folder + "Heat100.dat";
		String elasticOutputFile = folder + "Deform100.dat";
		String XYLineFile = folder + "Line100.dat";
		//String stressFile = folder + "Stress100.dat";
		//String LineStressFile = folder + "StressYLine100.dat";
		String expansion = folder + "Expansion100.dat";
		String temp00File = folder + "Temp000_100.dat";

		Mesh2Q4 mesh2D = new Mesh2Q4().square2D(24, 6);
		mesh2D.scale(1.0, 0.25);
		
		//mesh2D.displayBoundarys();
		//System.exit(0);

		//二维标量场有限元空间
		FES2Q41 fs = new FES2Q41(mesh2D);
		SRMatrix LMatrix = new SRMatrix(fs.GetNdof());
		SRMatrix RMatrix = new SRMatrix(fs.GetNdof());
		fs.assembleHeatThetaScheme(new double[] { kxx, density * specificHeat, dt, theta }, 
				                  new BVPType[] { BVPType.AXISYMMETRIC, BVPType.AXIS_CONSIST_MASS},
				                  LMatrix, RMatrix);
		//System.out.println(LMatrix);
		double[] x = new double[fs.GetNdof()];
		double[] RHS = new double[fs.GetNdof()];

		//二维弹性力学有限元空间
		FES2Q42 fs2 = new FES2Q42(mesh2D);
		SRMatrix K = new SRMatrix(fs2.GetNdof());
		fs2.assembleStiff(new double[] { E, Nu }, BVPType.AXISYMMETRIC, K);
		
//		K.sort();
//		System.out.println(K);
//		System.exit(0);
		
		double[] d = new double[fs2.GetNdof()];
		double[] heatLoad = new double[fs2.GetNdof()];
		double[] totalLoad = new double[fs2.GetNdof()];

		double[][] uv = new double[2][fs.GetNdof()];

		mesh2D.toTecplot(heatOutputFile, 0.0, x);
		mesh2D.toTecplotDeform(elasticOutputFile, 0.0, uv, 0.0);
		Tecplot.lineXYTime(XYLineFile, lineX, lineY, 0.0);
		//Tecplot.lineXYTime(LineStressFile, lineX, lineY, 0.0);

		double time = 0.0;
		fs.applyBC_MBN(LMatrix, Direct.All, 2);
		
		NewIterSSolver solverHeat = new NewIterSSolver(LMatrix);
		//DirectSSolver solverHeat = new DirectSSolver(LMatrix, fs.getnBand());

		fs2.applyBC_MBN(K, Direct.Y, 1, 3);
		fs2.applyBC_MBN(K, Direct.X, 4);
		
		NewIterSSolver solverElastic = new NewIterSSolver(K);
		//DirectSSolver solverElastic = new DirectSSolver(K, fs2.getnBand());

		double[][] vecX = new double[1][];
		ScalarFunc betaFunc = (xy, llabel, param) -> beta;
		ScalarFunc[] heatLoadFunc = { betaFunc };
		FileWriter fw2 = null;
		FileWriter fw1 = null;
		try {
			fw1 = new FileWriter(expansion);
			fw2 = new FileWriter(temp00File);
			fw1.write("0.0\t0.0\n");
			fw2.write("0.0\t0.0\n");
		} catch (IOException e) {
			// TODO Auto-generated catch block
			e.printStackTrace();
		}
		while (time < endTime) {
			time += dt;
			System.out.println("==================================================");
			System.out.println("time: " + time);
			fs.applyBC_MBN(RHS, Direct.All, 1.0, 2);
			//System.out.println(Arrays.toString(RHS));
			//System.exit(0);
			solverHeat.PCGSSOR(RHS, x, 1.5, 1);
			//solverHeat.PLUSolve(RHS, x, 1);
			mesh2D.toTecplot(heatOutputFile, time, x);
			Arrays.setAll(lineY, i -> fs.getValue(x, new double[] { lineX[i], 0.0}));
			Tecplot.lineXYTime(XYLineFile, lineX, lineY, time);

			// call fs.valueXY(x, 0.0D0, 0.0D0, val)
			// write(40, "(2F30.10)") time, val

			RMatrix.mul(x, RHS);

			Arrays.fill(heatLoad, 0.0);
			vecX[0] = x;
			fs2.assembleHeatLoad(heatLoadFunc, null, vecX, BVPType.AXISYMMETRIC, heatLoad);
			//System.out.println(MVO.toString(heatLoad));
			//System.exit(0);
			fs2.applyBC_MBN(heatLoad, Direct.Y, 0.0, 1, 3);
			fs2.applyBC_MBN(heatLoad, Direct.X, 0.0, 4);
			
			//System.out.println(MVO.toString(heatLoad));
			//System.exit(0);
			//solverElastic.PLUSolve(heatLoad, d, 1);
			solverElastic.PCGSSOR(heatLoad, d, 1.5, 1);
			//System.exit(0);
//
			fs2.extractComponent(d, uv);
//
			mesh2D.toTecplotDeform(elasticOutputFile, time, uv, 1.0 / ((1 + Nu) * alpha));
			// call fs2.computeStrain(u, v, ex, ey, exy)
			// call fs2.computeStress(ex, ey, exy, u, E, Nu, tx, ty, txy, u, 1)
			// tz = Nu * (tx + ty)
			System.out.println("Temperature: ");
			System.out.println("T min: " + MVO.min(x) + "\t " + MVO.max(x));
			System.out.println("Displacement: ");
			System.out.println("u min: " + MVO.min(uv[0]) + "\t " + MVO.max(uv[0]));
			System.out.println("v min: " + MVO.min(uv[1]) + "\t " + MVO.max(uv[1]));
			try {
				fw1.write(time + "\t " + fs.getValue(uv[0], new double[] {1.0, 0.0}) * 1.0 / ((1 + Nu) * alpha) + "\n");
				fw2.write(time + "\t " + fs.getValue(x, new double[] {0.0, 0.0}) + "\n");
			} catch (IOException e) {
				e.printStackTrace();
			}
			
			long end = System.nanoTime();
			System.out.println("Elapse time: " + (end - begin) / 1.0e9 + " s");

		}

	    try {
	    	fw1.flush();
			fw1.close();
			fw2.flush();
			fw2.close();
		} catch (IOException e) {
			e.printStackTrace();
		}
	}

}
